arXiv:1506.01913vl [math.NA] 4Jun2015 


Energy stable interior penalty discontinuous Galerkin 
finite element method for Cahn-Hilliard equation 


Biilent Karasozen^, Ayse Sanaydm Filibelioglu^’*, Murat Uzunca^ 

^Department of Mathematies & Institute of Applied Mathematies, Middle East Teehnieal 

University, 06800 Ankara, Turkey 

^Institute of Applied Mathematies, Middle East Teehnieal University, 06800 Ankara, 

Turkey 


Abstract 

An energy stable conservative method is developed for the Cahn-Hilliard 
(CH) equation with the degenerate mobility. The CH equation is discretized 
in space with the mass conserving symmetric interior penalty discontinuous 
Galerkin (SIPG) method. The resulting semi-discrete nonlinear system of or¬ 
dinary differential equations are solved in time by the unconditionally energy 
stable average vector field (AVF) method. We prove that the AVF method 
preserves the energy decreasing property of the CH equation. Numerical re¬ 
sults confirm the theoretical convergence rates and the performance of the 
proposed approach. 
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1. Introduction 

In this work, we consider Cahn-Hilliard (CH) equation in mixed form as 
a system of second order differential equations in a bounded domain Q C 
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or with periodic boundary conditions, where A denotes the Laplace operator 
and ii{u) is the non-negative mobility function. The parameter e provides a 
measure of the width of inter-facial layer, capturing the dominating effect of 
the reaction kinetics and represents the effective diffusivity. 

The CH equation is the most known model for phase separation and 
coarsening phenomena in a melted alloy [ 3 ] . It was also used as diffuse inter¬ 
face model for problems from fluid dynamics, material science and biology 
jl]. The variable u denotes the concentration of one of the species of the 
alloy and it is known as the phase state between materials. 

The CH equation describes a gradient flow with energy dissipation 


Ut = 


seju) 

Su 


with respect to the H inner product 


( 2 ) 


{u,v)h-i := {u,A ^v). 

In (|^ , the term S6 (u) / 5u stands for the variational derivative of the Ginzburg- 
Landau energy functional S given by 

£{u) = + F{u)^ dn, (3) 


with the potential function F(u) satisfying F'(u) = f{u). In the literature 
two different types of potential functions are used: the convex double-well 
potential 


F{u) 


(1 — 
4 


( 4 ) 
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and the non-convex logarithmic potential PE] 

9 9 

F{u) = - [n Inn + (1 — n) ln(l — n)] — (5) 

with 0 < 0 < 9c, where 9c is the transition temperature. For temperatures 9 
close to 9c, the logarithmic potential functional Q is usually approximated 
by the quartic double-well potential Q . Other types of logarithmic potential 
functions can be found in mw- Both potential functions Q and ([^ satisfy 
the monotonicity and the Lipschitz continuity conditions j9| 

(/(“l) - /(«2 ))(Mi - U 2 ) > -Ci{Ui - U2f, 

\f{ui) - fiu2)\ < Lf \ui - U 2 I , 
l/'(“i) - /'(“2)| < Lf \Ui - U 2 I , 


for ui,U 2 e with the constraints |ui^ 2 | < 1 and with Ci,Lf,Lf/ > 0. The 
two important properties of the CH equation are the conservation of mass 
under Neumann/periodic boundary conditions pQEllll 

/ u{t)dx = / u^dx, 

JQ, Jn 

and the decrease of total energy monotonically as t ^ 00 

In short, the CH equation describes a conservative gradient flow. The non¬ 
negative mobility function fj,(u) can be constant or degenerate [TU], where 
the diffusion process is restricted to the interface zone, i.e. u G [—1,1] and 
it is zero outside. Commonly adapted versions of the degenerate mobility 
functions are /i(u) = f3u{l — u) P E] and ii{u) = /3(1 — P P Bj for a 
positive jd. 

From numerical point of view, it is desirable to develop numerical schemes 
which preserve the energy decreasing and mass conservative properties of the 
CH equation. Energy stability implies that the total energy of the fully dis¬ 
crete CH equation dissipates in time analogously to the continuous energy 
Q. The schemes that preserve the discrete versions of the continuous en¬ 
ergy lead to approximate solutions which behave qualitatively similar to the 
continuous ones. Explicit methods are not suitable for time discretization of 


jji{u) I dVL. 
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the CH equation because they are not energy stable and require very small 
time steps due to the stability restrictions. 

In the literature, the CH equation with the constant mobility function has 
been discretized in space using finite differences dn, finite elements [I] and 
spectral methods [T]. Energy stable time discretization methods are based 
either on the convex splitting of the energy functional [2] or by adding stabi¬ 
lization terms to the energy functional F[u) [I2]. On the other hand, the CH 
equation with degenerate mobility has been discretized by continuous finite 
elements PC!, by local discontinuous Galerkin method piaiHi, discontin¬ 
uous Galerkin method with elements and with mixed finite elements P, 
finite differences [T3| and by spectral methods [H]. 

In this work, we use the mass conserving symmetric interior penalty dis¬ 
continuous Galerkin (SIPG) finite element method [El [E| for the space 
discretization. Since the discontinuous Galerkin method uses piecewise poly¬ 
nomials, which are discontinuous at the interfaces, the discontinuous Galerkin 
(DG) approximation allows to capture the sharp gradients or singularities for 
small e that affect the numerical solution locally. For the time integrator, it is 
well known that the first order backward Euler method is energy stable, i.e., 
the discrete energy decreases without any restriction of the step size At for 
very stiff gradient systems with e —?• oirzi. The only second order implicit en¬ 
ergy stable method is the average vector field (AVF) method which preserves 
the energy decreasing property of the gradient systems and the systems with 
Lyapunov functions [EllIZI. We prove that the energy decreasing property 
of the fully discrete GH equation with the degenerate mobility is preserved 
using SIPG in space and AVF method in time. Numerical convergence rates 
and numerical experiments for two dimensional problems with double well 
0 and logarithmic ([^ potential functions demonstrate the performance of 
the SIPG space discretization and the AVF time integrator. 

The remainder of the paper is organized as follows. In the next sec¬ 
tion, Section the SIPG discretization in space for the CH equation ([I]) 
is described and the semi-discrete system of ODEs in matrix-vector form is 
introduced. The fully discrete system using the AVF method as a time inte¬ 
grator is given in Section and in Section the energy decreasing property 
of the fully discrete system is proven. Numerical results are presented in Sec¬ 
tion to demonstrate the accuracy of the numerical approach. The paper 
ends in Section with some concluding remarks. 


4 


2. Semi-discrete formulations 


In this section, we outline the semi-discretization of the CH equation ([I]) 
using symmetric interior penalty discontinuous Galerkin (SIPG), a type of 
discontinuous Galerkin (DG) methods, space discretization. To be being, the 
classical (continuous) weak formulation of the CH equation 0 reads as: for 
a.e. t G (0, T] 


{ut,v)a + {ii{u)Vw,Vv)a = 0, e (6a) 

{w, v)n - (/(«), u)n - e^(Vu, Vu)n = 0, Vv e (^)) (6b) 

(u(0),v)n = (uo,v)n, Vv e (6c) 

where (•, •)q denotes the usual L^-inner product over the domain G. It is 
well-known that under certain regularity assumptions, the system 0 has a 
unique solution in the space L^(0,T; H^(Q)). 

In the sequel, we briefly describe the SIPG space discretization, and then, 
we introduce the matrix-vector form of the semi-discrete system of the CH 
equation ([^ as a system of ODEs. 

S.l. Discontinuous Galerkin discretization 

In this work, we use as a space discretization method for the CH equation 
([^ the SIPG method [I5l[l6] which is a type of discontinuous Galerkin finite 
elements method. Different from the classical finite elements methods, DG 
methods are suitable in the use of non-conforming grids, and requires lower 
regularity assumptions. 

Let 7/i be a family of shape regular (triangular) elements {Ki\ G Th 
such that D = U/L,, Ki n Af,- = 0 for Ki, Kj & Th, i ^ j- The diameter 
of an element K and the length of an edge E are denoted by Kk and Jie, 
respectively. We set the space of discontinuous test and trial functions 

Vh = {ue L\9) : u\k e W{K) VK G T] <t (7) 

where F^(K) denotes the set of all polynomials on K £ Th of degree at most 
q. We split the set of all edges Eh into the set £1° of interior edges and 
the set E^ of boundary edges so that Eh = E^U E^. Let the edge E be a 
common edge for two elements K and £7®. Since the functions u £ Vh are 
discontinuous along the inter-element boundaries, there are two traces of u 
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along E, denoted by u\e from inside K and E'\e from inside K^. Then, the 
jump and average of u across an interior edge E are dehned, respectively, by 

[n] = uIeHk + {«} = 

where uk and n^e denote the outward unit normal vector to the boundary 
of the elements K and on the edge E, respectively. Similarly, for a vector 
function Vu, the jump and average across an interior edge E are given by 

[Vu] = S/uIe • Ur + • liRe, {Vm} = -(Vu|e + VE\e)- 

On a boundary edge E C dQ, we set {«} = u and [u] = un, where n is the 
outward unit normal vector to the boundary. 

Then, in space SIPG discretized semi-discrete formulation of the CH 
equation 0 with a variable mobility function /r := jji,{u) reads as: set 
Uft,(0), w/i(0) G 14 be the projections (orthogonal L^-projections) of the initial 
conditions onto 14, find Uh{t),Wh{t) G 14 such that for almost every 

t G (0, T] and for all u G 14, we have 


{dtUh,Vh)Q + tthifi; Wh, Vh) = 0, (8a) 

{wh, Vh)a - {fiuh),Vh)a = ah{e^; ur vh)-, (8b) 


where the bi-linear (in last two argument) form Ofe(K; w, v) is given by 


ah{K;w,v 



(9) 


where the parameter a is called the penalty parameter and it should be suf¬ 
ficiently large to ensure the stability of the SIPG discretization as described 
in [16] with a lower bound depending only on the polynomial degree [Hj. In 
numerical experiments we take a = 3q{q + 1). 


2.2. Semi-discrete system in matrix-vector form 

Let Uh{t) and Wh{t) be the time dependent solutions of the semi-discrete 
system (|^. Then, Uh{t) and Wh{t) are of the form 

N riq N riq 

Uh{t) = E. Mt) = 

m=l j=l m=l j=l 
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where are the basis functions spanning the space f4, and and 
are the unknown coefficients. The number Ug stands for the local dimension 
depending on the basis polynomial order q and N is the number of (triangu¬ 
lar) elements. By substituting the expansions (10) into the system (|^ and 
choosing the test functions v = i = 1,... ,nq, k = 1,... ,N, we obtain 
the semi-linear systems of ordinary differential equations 


— 0 ) 

A,C + b{0-MC = 0, 

with the ordered unknown coefficient vectors and the basis functions 



= = 


• ? Sriq 5 SI ? • • 

• 5 ^Uq 5 • • 

cN 

• )Si ) • • • 

.N)T 

^^UqJ 5 

c 

= = 

(C!,.. 

• 5 Suq 5 Si ? • • 

• ’ Sriq ’ • • 


■ 5 Sriq/ 5 

P : 

= = 

(a},-- 


* * "> ^Uq") * 


\T 

' ’ ’ ^riq) 


In 0 , M denotes the mass matrix with the entries Mjj = 1 < 

bi ^^g X N, An and A^ are the stiffness matrices with the entries = 

and (Ae)ij = ah{e^]ip> ,p^), 1 < i,j <ngXN, and b is the 
non-linear vector of unknown coefficient vector ^ with the entries bi{^) = 
{f{uh), l<i<ngX N. 


3. Fully discrete system 

In this section, we give the fully discrete formulations of the CH equation 
0 through the time integration of the semi-linear system of ODEs 0. 
To do this, we consider the uniform partition 0 = to < ti < ... < tj = T 
of the time interval [0,T] with the uniform time step-size At = tk — 4-i, 
k = l,2,...,J. 

In this work, as a time integrator we use the average vector field (AVF) 
method which is a structure preserving time integrator for gradient systems 
[18] . For a given system of ODEs y — —VUiy), the AVF method is given by 

Vn = Vn-l - At j VV {ryn + (1 - T)yn-l)dT. (12) 

Jo 

The AVF method possesses the energy decreasing property without restric¬ 
tion of the step sizes At and is second order accurate in time. It represents a 
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modification of the implicit mid-point rule and for quadratic potentials U (y ), 
the AVF method reduces to the mid-point rule. Higher order variants of the 
AVF methods for Hamiltonian and Poisson systems with Gauss-Legendre col¬ 
location points are given in nzi. As Gauss-Legendre Runge-Kutta methods, 
the AVF method and higher order versions do not have damping property for 
very stiff systems, whereas for discontinuous Galerkin-Petrov methods and 
Radau H Runge-Kutta methods, the energy decreases monotonically with¬ 
out restriction of the step size At and the Lipschitz constant for VU{y). 
However, they require the solution of coupled system of equations, which in¬ 
creases the computational cost for two and three dimensional GH equations, 
where efficient solution techniques are required (TH]. 

Let for t = 0, Uhfi, Wh,o G 14 be the projections (orthogonal L^-projections) 
of the initial conditions wq onto 14 , and let rjo = (^ 0 ;Co)^ be the corre¬ 
sponding initial coefficient vector each of components satisfying (10). At 
a specific time t = tn, we denote the coefficient vector of the approxi¬ 
mate solutions {uh,niWh,n)'^ by rin = (^n)Cn)^- Then, the application of 
the AVF method (12) to the semi-discrete system leads to solving for 
n = 0,l,...,J — 1 the non-linear system of equations 


1 - 

1 _ 


1 

0 


1 - 

1 

LP. -i"J 

_Cn+l_ 


Jo K^^n+l + (1 - r)^n)dT_ 




(13) 

where we used the fact that the AVF method reduces the mid-point rule 
for linear terms. Further, the solution dependent mobility function y{uh) in 
the bilinear form is computed explicitly as for continuous finite elements 
in [H], i.e., on each time interval {tn,tn+i\ the mobility function is taken as 
fi = yiuh^n) where Uh^n is the known approximate solution from the previous 


time step. The system (13) can be written as a residual equation R{rf) = 
{Ri{r)), R 2 {ri))'^ = 0 with 


RliVn+l) — M{^n+1 — ^n) + A^[C,n^i + C,n) 

R2{r]n+l) = ^Ae(Cn+l - ]^M{C,n+l + Cn) + j ?>(TCn+l + (1 “ r)Cn)dT. 

(14) 

We solve the nonlinear system of equations ( [Idl ) using the Newton’s method. 
Starting with an initial guess rj^+i ~ (^n+n Cn^)^) i'be k—th Newton iteration 
to solve the nonlinear system of equations ([l4|) for the unknown vector r]n+i = 
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(^n+i,Cn+i)^ reads as 


js<‘> = 


>7l.+l'’ = ’ln+1 + S*‘’ 


A: = 0,1,... 


(15) 


until a user defined tolerance is satisfied. In (15), s = (si,S 2 )^ is the incre¬ 
ment, and J stands for the Jacobian matrix 


f M — A \ 

■^ = ( 1 ^, + 4 -Im)’ ^ + <' - ^Kn)dT ( 16 ) 

where + (1 — 'r)^n) is the Jacobian of the nonlinear form b{^) w.r.t. 

^ at ^ = T^n+i + (1 — 'r)^n- At each Newton iteration, the integral term in 
Jb is approximated by the fourth order Gaussian quadrature formula. 


4. Energy stability 


It is expected that fully discrete energy stable schemes should preserve 
the discrete energy dissipation as their continuous parts, which leads to qual¬ 
itatively better approximations. The continuous (in time) energy of the semi¬ 
discrete CH equation is given by [20] 


fV) = ^ l|VM|li,|„ + (n**). l)n. (17) 

On the other hand, the discrete DG counterpart of the continuous energy 
0 at a time = nAt reads as 

4gK) = ^ + (W, l)n 

+ X] K])i? + 


ae 

2hf 


(K],K])i 




In this section, we show that the AVF method applied to the semi-discrete 
system is energy stable through the discrete energy ( [I^ . Applying the 
AVF method to the semi-discrete system (1^ , we obtain for any z/, d G 14 


[u 


yjn+l _|_ yjU 


n+1 


At 


^)a + ^«fe(/^(«"'); = 0, 




1 


a 


- / (/(ru”+' + (l-T)u’^),d)adr=^a^e2;u”+i + u4d). 

(19) 
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Taking u = + w'^)l2 and d = — n” in (19), we obtain 




y^n+l _|_ 


+ ^ah{ij{u^);w^+^+ w^,w^+^+ w^) = 0, 


/,,,n+l _|_ „ .n \ pi 

^-J (/(Tn”+i + (l-rK),n"+i-n")odr = 


2. ^,n+l _|_ ^,n+l _ 


By using the identity {a + b,a — b)^ = {a^ — 6^, 1)^ and the bi-linearity of an, 
we get 


- u^, 


^n+l _j_ yjU 


+ ^atWti”);a)“+‘ + + «i”) = 0, 

( 20 ) 


w 


Ti'+l 


, 1 - I {f{TU^+^ + (1 - r)n”),- u^hdr = 


n+1 


n\ n+1 n,'bi\ 


Q ^0 




( 21 ) 


Using the Taylor expansions of the terms F[u^) and F(n”+^), and neglecting 
the higher order terms, we get 

F{u^) F{tu^+^ + (1 - r)n”) - /(r«”+^ + (1 - t)u^){t{u^+^ - u^)), 

F(n”+i) F{tu^+^ + (1 - t)u^) + firu^^^ + (1 - t)u^){1 - r)(n”+^ - u^)). 


Subtracting F{u'^) from F{u^~^^) leads to 

F{u^+^) - F{u^) ^ /(rn”+^ + (1 - t)u^)(u^+^ - u^) 

(F(n-+i), l)o - (F(n-), l)o ~ (/(rn"+i + (1 - t ) u ^), u ^+^ - 


((F(n-+i),l)n-(F(n-),l)n)dT 


/\/0 

Jo 


+ (1 - r)n”),n”+' - n”)adr 


n\ n+1 n,'F\ 


(F(n-+0,l)n-(i^K),l)n~ / (/(T«’^+' + (l-rK),n”+'-«”Mr 

Jo 

( 22 ) 
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We note that the bi-linear form ah satisfies 


ah{fiiu^);w^+^ + w^, + w^) = /rK) || VK+^ + w^) 

- 2 /" + w^)}[w^+^ + w^]ds 

EeEl 

+ E ^^11 + ”’"111V)' (23) 

egeo ^ 


Since all the terms in (23) are non-negative (see m Sec. 2.7.1] for positivity 

w^+^-\-w^) > 0. Similarly, 


of edge integral term), we have ah{ta{u)] 


we have ah{e 




u 


n+1 


) > 0 , ah{e^]u^,u^) 


subtracting (20) from (|21|) and substituting 


> 0. Using this identities, 
we obtain 


At 

T 


a,j (//(«”) 


+ w , w 


n+1 


+ w^) 


(F(u”+^), l)n + ia/,(e^; 

((F(u-),l)n + ia,(e2;u-,u-)^ < 0, 


which implies that < 8^q{u^). Hence, the AVF discretized 

scheme is energy stable through the discrete energy (18). 


5. Numerical results 

In this section, we present a set of numerical examples to confirm the 
accuracy, energy stability and mass conservation of our approach. All ex¬ 
amples are considered in two dimensional spatial domain with constant and 
degenerate mobility function, double-well (j^ and logarithmic ([^ potential 
functions under Neumann and periodic boundary conditions. 

5.1. Constant mobility and double-well potential under Neumann boundary 
conditions 

We consider 2D CH equation with the constant mobility function fj,(u) = 
1, diffusivity e = 0.1, and double-well potential (Q under the homogenous 
Neumann boundary condition [1]. The The exact solution in D = [—1,1] x 
[—1,1] for t G [0,1] is given by 

u{x, y, t) = cos(7rx) cos(7ry). 
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We add a load vector to the system so that the exact solution above solves 
the system. The errors and the numerical order of accuracy at time t = 1 
using the time step At = 0.5Ax for linear and quadratic DG polynomials 
are given in Table Both DG approximations satisfy the (k + l)-th order 
of accuracy for elements. 


Table 1: Example |5.1| Accuracy test with constant mobility and d ouble-well potential. 

Ax Dof L^-Error Order 


1/2 

24 

3.347 

- 

pi 1/4 

96 

1.633 

1.04 

1/8 

384 

4.810E-01 

1.76 

1/16 

1536 

1.079E-01 

2.16 


1/2 

48 

6.694E-01 

- 

p2 1/4 

192 

2.685E-01 

1.32 

1/8 

768 

3.376E-02 

2.99 

1/16 

3072 

3.733E-03 

3.18 


5.2. Degenerate mobility and double-well potential under periodic boundary 
conditions 

The next example is the CH equation in [2] with the exact solution 
m(x, y, t) = sin(x) sm{y) 

in the domain fl = [0, 27r] x [0, 27r] for t G [0,1] under periodic boundary 
conditions with the degenerate mobility function fj,{u) = 1 — and double¬ 
well potential 0- The effective diffusivity is taken as e = 1. For first 
order DG polynomials we use time step At = 0.003271 and for quadratic 
DG polynomials At = 0.000327r. Order reduction is observed in Table for 
linear and quadratic DG polynomials which is due to the nonlinearity of the 
degenerate mobility function. 

5.3. Constant mobility and double-well potential under Neumann boundary 
conditions: spinodal decomposition and nucleation 

We consider CH equation with constant mobility fi(u) = 1 and double¬ 
well potential function Q under homogenous Neumann boundary condition 
j2]. The computational domain is taken as D = [0,1] x [0,1] with e = 10“^. 
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Table 2: Example |5.2[ Accuracy test with degenerate mobility and dou ble-well potential. 

Mesh Size Dof L?' error Order 



t: 12 

24 

2.054 

- 

pi 

7r/4 

96 

5.742E-01 

1.84 


tt/S 

384 

1.566E-01 

1.87 


7r/16 

1536 

5.478E-02 

1.52 


7r/2 

48 

4.342E-01 

- 

p2 

7r/4 

192 

1.136E-01 

1.93 


tt/S 

768 

1.713E-02 

2.73 


7r/16 

3072 

4.895E-03 

1.81 


This problem represents the two main separation mechanisms: spinodal 
decomposition and nucleation. Both mechanisms in the CH equation are 
defined by the initial condition Uq{x) = r where h is a constant and r is 
a random number uniformly distributed on [0.005,0.005]. 

For h = 0, spinodal decomposition will be the governing separation mech¬ 
anism (see Fig. [^. The mixture separates from a randomly perturbed homo¬ 
geneous state [u = 0) giving rise to a striped pattern of complicated topology 
that coarsens over time. If we let the simulation evolve, the stationary so¬ 
lution would be a fully separated flow with two rectangular patches. While 
h 7 ^ 0, the separation mechanism is nucleation (see Fig. [^. We show the 
solution for h = 0.4. In the nucleation mechanism, isolated nuclei come up 
from the mixture. Again, the spatial micro-structure of the mixture coarsens 
over time. 

In both cases the discrete energy dissipates and the mass is conserved. 
Our results are similar to those in |2T], where for spatial discretization local 
discontinuous Galerkin method and for time discretization implicit convex 
splitting have been used. 

5.4- Degenerate mobility and logarithmic potential under Neumann boundary 
conditions 

We consider 2D CH equation with degenerate mobility function fi{u) = 
u{l — u) and the logarithmic potential function 

F{u) = 3000(n In n -|- (1 — n) ln(l — n)) -|- 9000m(1 — u) 

under homogenous Neumann boundary conditions [2] with the diffusion con¬ 
stant e = 1 in the domain D = [—0.5,0.5] x [—0.5,0.5] for t G [0,0.2]. The 
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initial condition is a random variation of uniform state u = 0.63 with a 
change no larger than 0.05. Fig. shows the evolution of the concentration 
held. The two phases in the concentration evolution, the phase separation 
stage and the coarsening process stage can be seen clearly. 

6. Conclusions 

We have presented numerical results for the CH equation under periodic 
and homogenous Neumann boundary conditions using SIPG discretization 
in space and AVF method in time. The numerical energy is decreasing in 
all examples and the numerical results are in good agreement to the those in 
the literature. 
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Figure 1: Example |5.3[ (Top-Middle) Solutions with constant mobility and double-well 
potential (spinodal decomposition, = 0), (Bottom) energy and mass evolutions. Linear 
DG elements and At = 10“^. 
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(e) 


(f) 


Figure 2: Example |5.3[ (Top-Middle) Solutions with constant mobility and double-well 
potential (nucleation, u ^0), (Bottom) energy and mass evolutions. Linear DG elements 
and At = 10“®. 
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(a) t = 0.000002 



X 


(b) t = 0.000008 



(c) t = 0.000128 


(d) t = 0.001 




Figure 3: Example |5.4| (Top-Middle) Solutions with degenerate mobility and logarithmic 
potential, (Bottom) energy and mass evolutions. Cubic DG elements and At = 10“^. 
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